#include "stdafx.h"
#include <iostream>
#include <stdlib.h>
#include <math.h>
#include <string.h> // memcpy
#include <fstream>

#include "calculationKN.h"

#ifndef M_PI
#define M_PI       3.14159265358979323846
#endif

#define FTX_VARIANT 1

double sigma = 1.0;
double k = 1.0;

double deltaX, deltaT;


double f_begin(double x){ // начальное условие, функция u(0, x)
	if (x <= 5)
		return sin(M_PI*x/5);
	else
		return 0;
}

double f_left(double time){ // левое граничное условие
	return 0;
}

double f_right(double time){ // правое граничное условие
	return 0;
}

double f(double t, double x){ // функция, задающая внешнее воздействие
	if ((x >= 5.0) && (x <= 5.2)) return 2;

	if (FTX_VARIANT == 1)
		return 0;
	else 
		if ((x >= 8.0) && (x <= 8.2)) return 2; else return 0;
}

void progonka(double** u, int count, int timeStep);

double **calculation(	
	int Nx, 	//количество интервалов по x
	int NTimes,     //количество интервалов по t
	double XEnd, 
	double TEnd
	)
{
	double XStart = 0.0;
	double TStart = 0.0;
	deltaX = (XEnd - XStart) / Nx;
	deltaT = (TEnd - TStart) / NTimes;

	double **data;
	int stepsX = Nx + 1;
	int stepsT = NTimes + 1;

	data = new double*[stepsT](); // ряды в рассматриваемые моменты времени
	for(int i = 0; i < stepsT; i++){
		data[i] = new double [stepsX](); // массивы значений в момент времени
	}

	double **u;
	u = new double*[2]; // на смену 2 временных ряда
	u[0] = new double [stepsX]();
	u[1] = new double [stepsX]();

	// считаем в начальное время:
	for (int i=0; i<stepsX; i++)
	{
		double x = XStart + i * deltaX;
		u[0][i] = f_begin(x);
	}

	memcpy(data[0], u[0], sizeof(double) * stepsX);

	for(int i = 1; i < stepsT; i++)
	{
		// строим электронные плотности во все моменты времени
		int curr = i % 2;	// вычисление в n+1 момент времени
		int prev = (i-1) % 2; // по n-ному моменту времени

		u[curr][0] = f_left(i * deltaT);
		u[curr][stepsX-1] = f_right(i * deltaT);

		progonka(u, stepsX, i);

		memcpy(data[i], u[curr], sizeof(double) * stepsX);
	}

	return data;
}

void progonka(double** u, int count, int timeStep)
{
	int curr = timeStep % 2;
	int prev = (timeStep - 1) % 2;
	//применяем метод прогонки к неизвестным u[curr][1] ... u[curr][stepX - 2] -> table[0] ... table[n-1]
	int n = count - 2;	//количество неизвестных
	double* table[2];
	table[0] = &(u[0][1]);
	table[1] = &(u[1][1]);
	//формируем столбцы A, B, C 3-х диагональной матрицы
	//и столбец F правой части
	//см. http://ru.wikipedia.org/wiki/%CC%E5%F2%EE%E4_%EF%F0%EE%E3%EE%ED%EA%E8

	//первый и последний элементы по-особому, т.к. там 2 неизвестных
	//см. http://ru.wikipedia.org/wiki/Трёхдиагональная_матрица
	double* A = new double[n];
	double* C = new double[n];
	double* B = new double[n];

	double* F = new double[n];

	//высчитываем коэффициенты A, B, C, F
	//теория 1 http://ipicture.ru/Gallery/Viewfull/33199189.html
	//теория 2 http://ipicture.ru/Gallery/Viewfull/33199191.html
	double r = sigma / (2*deltaX*deltaX); // вспомогательный коэффициент
	A[0] = 0;
	C[0] = (1 / deltaT + 2*r + k);
	B[0] = -r;
	F[0] = r*(table[prev][1]-2*table[prev][0]+table[prev][-1]+table[curr][-1])+table[prev][0]*(1/deltaT) + f( (timeStep-1) * deltaT, 1*deltaX);
	//
	A[n-1] = -r;
	C[n-1] = (1/deltaT + 2*r + k);
	B[n-1] = 0;
	F[n-1] = r*(table[prev][n]-2*table[prev][n-1]+table[prev][n-2]+table[curr][n])+table[prev][n-1]*(1/deltaT)+f( (timeStep-1) * deltaT, n * deltaX);

	for (int i=1; i<n-1; i++)
	{
		A[i] = -r;
		C[i] = (1/deltaT+2*r+k);
		B[i] = -r;
		F[i] = r*(table[prev][i+1]-2*table[prev][i]+table[prev][i-1]) + table[prev][i]*(1/deltaT)+f( (timeStep-1) * deltaT, (i+1) * deltaX);
	}

	//doProgonka(A,C,B,F)
	{
		double* alpha = new double[n];
		double* beta = new double[n];

		alpha[0] = beta[0] = 0;
		for (int i=0; i<n-1; i++)
		{
			alpha[i+1] = -B[i] / (A[i]*alpha[0] + C[i]);
			beta[i+1] = (F[i] - A[i]*beta[i]) / (A[i]*alpha[i] + C[i]);
		}

		table[curr][n-1] = (F[n-1] - A[n-1]*beta[n-1]) / (C[n-1] + A[n-1]*alpha[n-1]);
		for (int i=n-2; i>=0; i--)
			table[curr][i] = alpha[i+1]*table[curr][i+1] + beta[i+1];

		delete []alpha;
		alpha = NULL;
		delete []beta;
		beta = NULL;
	}
	//
	delete []A;
	A = NULL;
	delete []B;
	B = NULL;
	delete []C;
	C = NULL;
	delete []F;
	F = NULL;
}

int f_main(int argc, char **argv)
{
	double XStart = 0.0;
	double XEnd = atof(argv[1]); // конец отрезка по x
	int Nx = atoi(argv[2]); // количество интервалов, например 100
	deltaX = (XEnd - XStart) / Nx; 

	double TStart = 0.0;
	double TEnd = atof(argv[3]); // конец времени моделирования
	int NTimes = atoi(argv[4]); // количество интервалов по времени
	deltaT = (TEnd - TStart) / NTimes;

	// double **data = calculation(deltaX, time, deltaTShow);
	double** data = calculation(Nx, NTimes, XEnd, TEnd);

	int stepsT = NTimes + 1;
	int stepsX = Nx + 1;

	std::cout << "{\"u\":[";
	for(int i = 0; i < stepsT; i++)
	{
		std::cout << "[";
		for(int j = 0; j < stepsX; j++)
		{
			std::cout << data[i][j];
			if(j < (stepsX - 1))
				std::cout << ",";
		}
		std::cout << "]";
		if(i < NTimes)
			std::cout << ",";
	}
	std::cout << "],\"f\":[";
	for(int i = 0; i < stepsT; i++)
	{
		std::cout << "[";
		for(int j = 0; j < stepsX; j++)
		{
			std::cout << f(i * deltaT, j * deltaX);
			if(j < (stepsX - 1))
				std::cout << ",";
		}
		std::cout << "]";
		if(i < NTimes)
			std::cout << ",";
	}
	std::cout << "]}";
	return 0;
}

double** run(int argc, char **argv)
{
	double XStart = 0.0;
	double XEnd = atof(argv[1]); // конец отрезка по x
	int Nx = atoi(argv[2]); // количество интервалов, например 100
	deltaX = (XEnd - XStart) / Nx; 

	double TStart = 0.0;
	double TEnd = atof(argv[3]); // конец времени моделирования
	int NTimes = atoi(argv[4]); // количество интервалов по времени
	deltaT = (TEnd - TStart) / NTimes;

	// double **data = calculation(deltaX, time, deltaTShow);
	double** data = calculation(Nx, NTimes, XEnd, TEnd);

	int stepsT = NTimes + 1;
	int stepsX = Nx + 1;

	freopen("calculation.txt","w",stdout);


	std::cout << "{\"u\":[";
	for(int i = 0; i < stepsT; i++)
	{
		std::cout << "[";
		for(int j = 0; j < stepsX; j++)
		{
			std::cout << data[i][j];
			if(j < (stepsX - 1))
				std::cout << ",";
		}
		std::cout << "]";
		if(i < NTimes)
			std::cout << ",";
	}
	std::cout << "],\"f\":[";
	for(int i = 0; i < stepsT; i++)
	{
		std::cout << "[";
		for(int j = 0; j < stepsX; j++)
		{
			std::cout << f(i * deltaT, j * deltaX);
			if(j < (stepsX - 1))
				std::cout << ",";
		}
		std::cout << "]";
		if(i < NTimes)
			std::cout << ",";
	}
	std::cout << "]}";
	
	fclose(stdout);

	return data;
};